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Preface 


The  purpose  of  this  study  was  twofold:  to  determine 
the  validity  of  two  new  silver  potentials  and  to  verify  a 
new  global  minimum  configuration  technique  for  silicon. 
First,  molecular  statics  was  used  to  test  two  new  Dagen 
silver  potentials.  In  particular  this  was  accomplished  by 
calculating  the  formation  energy  of  vacancy  point  defects  in 
the  face-centered  cubic  metal  silver  via  computer 
simulation.  Then,  by  comparing  the  results  obtained  with 
known  experimental  data  the  viability  of  the  new  potentials 
versus  the  prior  version  was  determined.  Second,  molecular 
statics  in  conjunction  with  molecular  dynamics  was  used  to 
determine  the  global  minimum  configuration  of  a  two 
dimensional  36  atom  cluster  using  a  modified  version  of 
Stillinger  and  Weber’s  silicon  potential.  The  resulting 
computational  times  were  then  compared  to  those  obtained  by 
utilizing  a  prior  method  propossed  by  Biswas  and  Haaman 
utilizing  Langevin  molecular  dynamics  techniques  (Biswas  and 
Haaman ,  1986 ) . 

In  performing  the  computer  simulation  and  writing 
this  thesis  I  thankfully  have  had  a  great  deal  of  support 
and  help  from  others.  I  am  deeply  indebted  to  my  faculty 
advisor,  Capt  M.  Sabochick,  for  his  continuing  patience  and 
ass istance . 
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Debra  L.  Richlin 
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Abstract 


Two  new  silver  potentials  derived  by  Dagens  were 
validated  and  a  new  annealing  technique  was  investigated  in 
this  thesis.  It  was  found  via  a  molecular  statics  technique 
that  the  two  new  silver  potentials  were  Jess  precise  than 
the  older  silver  potential.  The  calculated  vacancy  of 
formation  energy  of  the  older  silver  potential  deviated  from 
the  experimental  value  by  21%  whereas  the  two  new  silver 
potential  values  differed  by  48%  and  64%  respectively. 

The  new  annealing  technique  using  a  combination  of 
molecular  statics  and  molecular  dynamics  was  found  to  be 
approximately  ten  times  faster  than  a  previous  method 
proposed  by  Biswas  and  Haaman  using  Langevin  dynamics. 

In  particular,  both  techniques  using  a  two-dimensional 
cluster  of  36  atoms  and  a  modified  version  of  the  Stillinger 
and  Weber  silicon  potential  were  tested  on  the  CRAY  computer 
system.  The  new  annealing  technique  acheived  a  local  minimum 
configuration  of  -45.3  eV  in  36  seconds  of  CRAY  CPU  time 
while  the  BH  Langevin  technique  took  over  360  seconds  of 
CRAY  CPU  time  to  find  the  same  local  minimum. 
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ATOMISTIC  SIMULATION  OF 
THE  VERIFICATION  OF  DAGENS ’  POTENTIALS  AND 
THE  ANNEALING  OF  ATOMIC  CLUSTERS 

I .  INTRODUCTION 

Atomistic  simulation  is  the  computer  simulation  of  a 
material  as  a  collection  of  individual  atoms.  Two  major 
types  of  atomistic  simulation  are  molecular  dynamics,  MD , 
and  molecular  statics,  MS.  Molecular  dynamics  is  used  to 
simulate  a  material  at  finite  temperatures.  Many  bulk 
properties  of  the  material  can  thus  be  calculated  such  as 
pressure,  internal  energy  and  enthalpy.  Molecular  statics 
is  used  to  determine  stable  and  metastable  atom  configura¬ 
tions  and  is  strictly  valid  only  at  0°  Kelvin. 

MD  and  MS  were  employed  in  completing  the  two  goals  of 
this  work.  These  goals  were  first  to  validate  Dagens’  new 
silver  potentials  and  second  to  investigate  a  possible 
faster  technique  for  finding  global  energy  configurations. 

In  atomistic  computer  simulations  the  material  in  ques¬ 
tion  is  modeled  by  individual  atoms  connected  via  inter¬ 
atomic  potentials.  These  interatomic  potentials  dictate  how- 
far  apart  the  atoms  can  be  and  still  react  with  each  other. 
They  also  dictate  the  strength  and  direction,  repulsive  or 
attractive,  of  the  interactive  forces  (Glasstone,  1979:253- 
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257).  Thus  to  simulate  a  real-world  material  a  reliable 
interatomic  potential  that  correlates  well  with  experimental 
data  is  necessarily  the  cornerstone  of  a  working  code. 

The  potentials  used  to  simulate  silver  were  three  noble 
metal  potentials  derived  from  first  principles  by  Dagens 
(Lam,  1983:2510).  The  first  and  oldest  of  these  potentials 
was  tested  by  Lam  (Lam,  1983)  to  verify  its  accuracy  with 
experimental  data.  It  was  found  that  for  copper  the 
potential  agreed  quite  well,  within  10%  of  the  experimental 
data  value  for  the  vacancy  energy  of  formation  (the  amount 
of  energy  required  to  form  a  vacancy  in  a  perfect  crystal). 
But  for  silver  the  energy  of  formation  was  over  estimated  by 
23%  when  compared  with  experimental  data  (Lam,  1983:2514). 
Subsequently  Dagens  revised  this  potential  in  an  attempt  to 
achieve  better  correlation  with  experimental  data.  The  two 
new  potentials  that  he  derived  were  tested  as  part  of  this 
project  using  a  molecular  statics  approach  to  find  the 
energy  of  vacancy  formations.  This  value  when  compared  with 
experimental  data  reflected  the  validity  of  the  two  new 
potentials . 

Biswas  and  Haaman ,  BH,  ( BH ,  1986:897-898)  utilized 

Lagevin  Molecular  dynamics  to  find  the  global  minimum  energy 
configuration  of  a  32  atom  cluster.  That  is  the  stable  or 
lowest  energy  configuration  of  the  system.  They  used  a 
technique  where  the  atoms  were  thermally  excited  to  the  melt 
point  and  then  allowed  to  cool  through  a  series  of  cooling 
stages  to  a  global  minimum  configuration.  At  each  stage  of 
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the  cooling  process  they  recalculated  the  positions  of  the 
atoms  using  Langevin  dynamics  in  conjunction  with  the  Stil- 
linger  and  Weber  silicon  potential.  To  achieve  their 
results  BH  had  to  calculate  the  Langevin  equation  200-300 
times  for  each  temperature  decrement  ( BH ,  1986:897-898). 

This  method,  although  faster  than  other  techniques,  is  still 
somewhat  computationally  time  consuming.  The  second  part  of 
this  paper  proposed  using  a  combination  of  molecular  statics 
and  molecular  dynamics  to  find  the  global  minimum  configura¬ 
tion  within  less  timesteps  than  those  required  using  the  3H 
Langevin  technique. 

In  their  article  BH  refer  to  their  technique  as  a  means 
of  simulating  the  annealing  of  a  cluster  of  atoms  (BH,1986). 
Experimentally  annealing  is  the  gradual  removing  of  energy 
from  a  heated  system  until  it  reaches  a  stable  configura¬ 
tion  ( Chalmers , 1982 :  109-  1  1 1 )  .  Through  a  slow  annealing 
process  many  defects  disappear  from  the  crystal  as  they  are 
given  time  to  recombine,  an  interstitial  with  a  vacancy,  or 
reposition  themselves  according  to  the  potentials  influence 
'Corbett,  1986:36).  For  the  purpose  of  this  paper,  the 
cyclic  cooling  process  used  in  the  computer  simulations  will 
be  referred  to  as  annealing. 

Overview . 

Chapter  2  describes  the  interatomic  potentials  used  in 
this  work.  The  concepts  of  two  body  versus  three  body 
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II .  INTERATOMIC  POTENTIALS 


Interatomic  potentials  are  used  to  determine  the 
interatomic  forces  and  energies  between  atoms.  These  forces 
and  energies  are  used  by  molecular  dynamics  and  molecular 
statics  to  calculate  both  micro-  and  macro-  properties  of 
materials  such  as  the  vacancy  energy  of  formation  and  the 
global  energy  configuration  of  an  annealed  solid.  The 
results  obtained  when  simulating  a  particular  material 
relies  heavily  on  the  accuracy  of  the  interatomic  potentials 
ability  to  effectively  model  the  real-world  particle 
interactions  (Sabochick,  1985:20). 

Two  and  Three  Body  Potentials. 

Generally  potentials  are  based  on  the  sum  of  one  body, 
two  body,  three  body  etc.  interactions  between  atoms.  (SW, 
1985)  Therefore  interatomic  potentials  can  be  classified 
according  to  the  number  of  atoms  that  interact  with  each 
other  at  any  one  time.  Materials  that  can  be  adequately 
portrayed  by  the  simultaneous  interaction  of  two  atoms 
utilize  the  two  body  term.  Materials  that  are  more  complex, 
whose  atoms  interact  with  each  other  in  triplets,  utilize 
both  the  two  and  three  body  terms.  The  latter  set  are 
referred  to  as  three  body  potentials  as  they  include  the 
three  body  term.  These  two  and  three  body  terms  can  be 
expressed  via  a  generic  formula. 


The  general  expression  is  of  the  form 


0 


=  £ 


Vi(i)  + 


£v*(i.j>  +  £v 


a  (  i  ,  j  ,  k  )  + 


(  1  ) 


i  <  i 


i  j  k 
i  <  >  <  k 


The  first  term  describes  external  and  wall  forces  that 
interact  with  the  atoms  (SW,  1983).  As  the  crystal  is 
normally  taken  to  be  isolated  from  external  forces  the  one 
body  term  is  usually  not  included  in  the  potential 
(SW,  1983). 

The  second  term  describes  the  interactions  between 
pairs  of  atoms  which  are  not  affected  by  the  presence  of 
other  neighbors.  This  two  body  term  is  commonly  known  as  a 
pair  potential.  Pair  potentials  are  the  most  widely  used 
potentials.  The  two  body  term  sufficiently  approximates 
some  materials,  is  easily  implemented  and  computationally 
faster  than  more  complex  potentials  (Thee,  1986). 

The  third  or  three  body  term  is  used  when  the  material 
has  strong  directional  bonds  (SW,  1983).  This  three  body 
potential  is  used  to  stabilize  the  system  and  provides  the 
unique  geometrical  configuration  that  is  an  integral  part  of 
certain  materials  such  as  silicon.  In  particular,  the  three 
body  term  forces  the  bond  angles  to  a  given  value.  This  term 
portrays  the  simultaneous  interaction  of  three  atoms  bonded 
together  by  their  individual  atomic  forces.  Implicit  in  the 
three  body  term  is  thus  an  angular  dependence  (Fig  1.) 
(Sabochick,  1985). 
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Figure  1.  Three  body  atomic  interation.  Dotted  lines 
represent  interatomic  force  fields  binding  the  atoms 
together.  The  o  represents  the  individual  atoms. 


Empirical  veraua  Non-empirical  Potentials. 

Interatomic  potentials  can  also  be  classified  according 
to  the  way  they  are  derived.  Two  classes  of  potentials  are 
empirical  and  non-empirical.  Empirical  potentials  are 
derived  by  fitting  a  specific  form  of  a  potential  such  that 
it  will  reproduce  experimental  data.  The  non-empirical 
potentials  are  derived  from  first  principles  employing  a 
quantum  mechanical  treatment  of  atomic  interactions 
(Sabochick,  1985:20-27). 

Although  empirical  potentials  reproduce  the  data  to 
which  they  are  fitted,  they  may  be  inaccurate  when  modeling 
other  parameters  than  the  ones  used  for  the  initial  data 
base.  An  example  of  an  empirical  three-body  potential  is 
the  Stillinger  and  Weber  potential  for  silicon  ( SW , 

1983)  . 
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One  type  of  non-empir ical  or  first  principles  potential 
is  derived  using  the  psuedopotential  method  (Sabochick,  1985 
:27-36).  These  potentials  have  a  number  of  common 
properties.  One  of  these  common  properties  is  that  beyond 
the  first  lattice  constant  the  potentials  demonstrate 
oscillatory  motion.  At  large  distances  these  oscillations 
usually  converge  to  periodic  Friedel  oscillations  (Lam, 1983, 
Sabochick , 1986  )  . 

The  Three  Interatomic  Pattens*  Potentials. 

Dagens ’  potentials  were  all  derived  from  first 
principles  using  the  psuedopotential  method  (Lam, 

1983:2504).  Figures  2,  3  and  4  portray  the  three  Dagens 
potentials  over  different  ranges  of  interatomic  distances. 
Figures  2  and  3  show  that  the  potentials  are  quite  disimilar 
at  r  <  a,  where  a  is  one  lattice  constant.  The  lattice 
constant  for  silver  is  4.087  Angstroms.  Note  in  Fig.  3  that 
the  second  new  potential  in  this  region  has  a  prominent 
negative  well.  The  old  potential  is  the  most  strongly 
repulsive  of  the  three.  The  first  new  potential  varies  only 
slightly  from  the  old  potential  in  this  region.  However 
viewing  Fig.  4  one  can  see  that  at  distances  greater  than  a 
few  lattice  constants  the  new  potentials  resemble  each  other 
while  the  old  potential  dominates  with  pronounced  Friediel 
oscillations.  The  Friedel  oscillation  dependence  which  is 
prevalent  in  psuedopotentials  can  thus  be  seen  in  Fig.  4. 
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Figure  3.  Datjens'  Three  Silver  Potentials  t  r  <  2a  )  . 
Where  a  is  the  lattice  constant  for  silver. 
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The  older  version  of  the  interatomic  Dagens  potential  is 


C  j  cos ( x  )  S  i  s  in ( x  ) 

0  =  {  1  -exp  [  -a  (  r  -  ri*)]J[(Ci+ - )  -  +  -  .  .  —  ] 

r 2  r 3  r 5 


+ 


n 


Baexp( -br  3  ) 


r  3  -  n 


+  D  exp(-gr) 


(2) 


where  x  =  Y  +  2  kfTf  (k  is  the  zero  order  Fermi  momentum 
for  noble  metals)  (Lam,  1983).  The  potential  0  is  in  atomic 
units  with  the  calculational  parameters  presented  in  Table 
I.  The  first  terra  in  the  above  equation  represents  Friedel 
oscillations  at  large  distances.  The  second  term  represents 
a  short  range  repulsive  force.  The  third  term  is  related  to 
the  forces  resulting  from  scattering  between  two  resonant  d 
states  (Lam, 1983).  The  old  Dagens  parameters  are  presented 
in  Table  I . 


TABLE  I 

Interatomic  Potential  Parameters  for 
Silver  (Old  Potential) 


a 

0.1 

b 

0. 109 

r  i 

3. 

Bo 

17.669 

Co 

0.265 

Bt 

-0.0459 

Ci 

-0.265 

Bj 

-0.2044 

s, 

-20.85 

B, 

-0.02728 

Y 

1 .429 

g 

0.8 

kr 

0.7189 

D 

0.2819 
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The  only  difference  between  the  two  new  Dagens 
potentials  are  the  values  of  the  constant  parameters.  Thus 
the  two  new  potentials  are  of  the  form. 

Ft  o  t  =  Fsp  +  Fdd  (3) 

where 

Fsp  (r)  =  (f ,  +  f 2  )  C  +  f a  (  1  -  C)  (4) 

where 


ft  =  At  exp( -B i  r)  cos(Kt  r  +  FI,  )  (5) 

C  =  l/(  1  +  exp(  2 ( r-c )  )  )  (6) 


and 


2 

Fdd(r)  =  [(DO  +  D2  r  )  exp(-GO  r)  +  D1  exp(-Gl  r ) I  C' 

+  D3  (1  -  C’)  cos  (  K3  (  r  -  RO  )  ) 

*  exp  (-G(r-RO))  (7) 

where 

C*  =  1/(1+  exp (  2 ( r-c  ’  )  )  )  (8) 

The  above  terras  involve  the  corresponding  resonant  s, 
p,  and  d  states  of  the  atom.  The  first  term,  Fsp,  is  a 
result  of  electronic  transitions  between  the  resonant  s 


state  at  one  ionic  site  and  the  resonant  p  states  at  a 
different  ionic  site.  The  second  term,  Fdd,  involves  the 


contributions  associated  with  the  scattering  between  two 


resonant  d  states.  Table  II  presents  the  common  parameters 
between  the  two  new  interatomic  Dagens  potentials,  while 
Tables  III  and  IV  present  respectively  the  FSp  portion  of 
Dagens  potentials  1  and  2  (Lam, 1986). 

TABLE  II 

Interatomic  Potential  Parameters  for 
Silver  (  F  part  of  New  Dagens’  Potentials) 
dd 


DO 

0.667611 

G 1 

1.763201 

D 1 

16.094525 

K3 

1.395939 

D2 

0.308763 

G3 

0.500000 

D3 

-0.000011 

c* 

7.238867 

GO 

1.481220 

R0 

10.427867 

TABLE  III 

Interatomic  Potential  Parameters  for 


Silver 

(  F  part  of 

sp 

New  Dagens 

Potent ial  1  ) 

A  i 

27.846879 

B  i 

1 . 301653 

A  i 

0.538025 

B, 

0.666321 

A  3 

-0.008047 

B, 

0. 405860 

k  i 

0.654843 

FI  i 

2.019130 

k  2 

-0.298960 

FI  2 

1  .745467 

k  3 

l  .  1745467 

FI  3 

-0.020818 

c 

7.9488 
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TABLE  IV 


Interatomic  Potential  Parameters  for 
Silver  <  F  part  of  New  Dagens *  Potential  2) 
sp 


A  1 

-24.090029 

B  1 

1 .249204 

A  2 

0.516413 

Ba 

0.649002 

A  3 

-0.010724 

b3 

0.460805 

k  1 

0.686321 

FI  1 

1 .932997 

k  2 

-0.029370 

FI  2 

0. 165145 

k  3 

1 . 177000 

FI  3 

0 .279729 

c 

7.9913 

The  Modified  Stillinger  and  Weber  Silicon  Potential. 

The  Stillinger  and  Weber  empirical  potential,  0,  is 
composed  of  a  two-body  term,  fz,  and  a  three-body  term,  f3, 
as  follows: 

r  i  i  r  i  r  j  rk 

0  =  ef  2  ( - )  +  ef  j  (  - ,  -  ,  -  )  (9) 

s  s  s  s 

where  r, j  is  the  radial  distance  between  atom  i  and  atom  j. 
The  terms  ri,  r3,  and  ru  are  respectively  atoms  i,  j,  and 
k’s  radial  positions,  e  and  s  are  energy  and  length  units. 
The  e  was  chosen  as  50  kcal/mole  to  give  fa  depth  of  -1  and 
s  was  chosen  as  0.20951  nm  to  make  f 2 ( 2  1  '  6 )  vanish.  Thus  a 
minimum  occurs  at  r  =  2 1  ' 6  and  the  largest  two  atom  binding 
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energy  is  -1.  ( Dodson , 1 986 )  The  f2  term  specifies  the 

energy  between  two  atoms  as  a  function  of  distance.  The  fj 
term  possesses  full  rotational  and  translational  symmetry 
(SW,  1983). 

The  Stillinger  and  Weber  two-body  term  is 


f  2  (  r  )  =  I 


A(B  r-P  -  r-<*)  exp[l/(r-a)  )  ,  r  <  a 

(  10  ) 

0  ,  r  >  a 


and  the  three-body  term  is  the  symmetrized  sum 

f  3  =  h(ijk)  +  h(jik)  +  h(jki)  (11) 

where  h(ijk)  depends  on  the  distances  r, ,,  r,k,  and  the 
angle  8 ; jk  (Dodson,  1983).  For  r  <  a  ,  h  has  the  following 
form 


h( r  i  j , r j k , 0 j  i  k  )  =  L'exp{Y[ l/(r ,  j  -  a)]  +  Y[ 1 / ( r  ,  u-a )  ]  ) 

*  [cos(8jik)  -  cos(8t)!  (12) 

The  parameters  that  Stillinger  and  Weber  used  in  the 
above  potentials  to  simulate  silicon  are  presented  in 
Table  V.  (SW,  1983).  With  only  a  few  exceptions  these  were 
used  as  shown  throughout  this  study.  Note  that  all  the 
values  are  given  in  reduced  units  where  s  =  2.0951  A°  and 
e  =  2.167  ev/atom  pair  ( Blai sten-Baro jas ,  1986). 

To  give  a  well  defined  hexagonal  pattern  in  two- 
dimensions  a  number  of  changes  in  parameters  were  necessary 
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TABLE  V 


Interatomic  Potential  Parameters  for 
Silicon  (  Stillinger  and  Weber  Potential) 


A 

7.049556277 

L’ 

21  . 

B 

0.602245584 

0  t 

-1/3 

q 

0.0 

P 

4  . 

a 

1 . 80 

Y 

1  .  20 

during  this  study.  One  change  shown  earlier  by  Dodson  was 
that  the  L’  term  needed  to  be  increased  from  21  to  a  larger 
value  (Dodson  chose  60)  to  adequately  model  silicon  crystal 
growth  ( Dodson , 1 983 ) .  The  L'  term  dictates  the  preference  in 
two-dimensions  of  the  hexagonal  structure  of  the  crystal 
over  other  atomic  configurations.  By  strengthening  the 
three-body  term  the  hexagonal  pattern  dominated  the  atomic 
cluster.  We  found  L’  =  30  to  be  adequate  for  the  purpose  of 

modeling  the  two-dimensional  melt  and  subsequent  annealing  A 

of  flie  simulation  material. 

Another  parameter  that  was  changed  during  this  study 
was  cos(0t)i  which  also  determine  the  resulting  geometry  of 
the  solid.  Stillinger  and  Weber  used  eos(0t)  =  -1/3  which 
gives  a  tetrahedral  geometry  in  three-dimensions.  We  chose 
to  use  a  120  degree  angle  which  will  give  a  hexagonal 
structure  in  two-dimensions.  The  parameters  used  during 
this  study  are  presented  in  Table  VI  with  the  two  modified 
parameters  changed  as  specified  above. 
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TABLE  VI 


Modified  Interatomic  Potential  Parameters 
for  Silicon  (  Stillinger  and  Weber  Potential) 


A 

7 .049556277 

L’ 

30  . 

B 

0 . 602245584 

9  « 

-0.5 

q 

o 

o 

P 

4  . 

a 

1 . 80 

Y 

1  .  20 
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ATOMISTIC  SIMULATION 


Atomistic  simulation  is  a  computer  technique  whereby 
a  material  is  simulated  via  atoms  connected  by  interatomic 
potentials.  These  potentials  are  used  to  provide  force  and 
energy  calculations  for  the  computational  methods  which 
adjust  the  atoms  positions.  These  computational  methods  can 
be  arranged  into  two  general  categories.  They  are  molecular 
dynamics  and  molecular  statics.  The  latter,  molecular 
statics,  was  used  to  verify  Dagens ’  silver  potentials.  Also 
a  combination  of  molecular  dynamics  and  molecular  statics 
was  used  to  investigate  a  new  annealing  technique.  A 
modified  version  of  molecular  dynamics  using  the  Langevin 
equations  was  also  used  in  an  annealing  technique  by  Biswas 
and  Haaman,  BH.  The  BH  technique  was  used  to  compare  the 
results  obtained  with  the  new  annealing  technique.  Both  of 
these  annealing  techniques  are  presented  in  the  latter  part 
of  this  chapter. 

The  following  paragraph  delineates  the  sections  covered 
in  this  chapter.  The  two  major  sections  covered  are  general 
computational  procedures  and  the  two  annealing  techniques. 
Under  general  computational  procedures,  molecular  dynamics, 
molecular  statics  and  a  modified  molecular  dynamics 
technique  using  Langevin  dynamics  is  covered.  Under  the  two 
annealing  techniques,  the  BH  method  using  Langevin  dynamics 
and  the  new  annealing  method  is  covered.  Each  annealing 
technique  is  outlined  along  with  particular  computational 
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tools  and  specific  parameters  utilized.  Finally  this 
chapter  concludes  with  a  description  of  the  constraint 
potential  used  in  the  BH  Langevin  and  the  new  annealing 
techn i ques . 


General  Computational  Procedures. 

Molecular  Dynamics. 

Molecular  dynamics  is  a  technique  whereby  one  solves 
Newton's  equations  of  motion  for  the  atom  positions  and 
velocities  of  a  system  at  finite  temperatures.  For  a  system 
of  M  atoms  with  velocities  Vi  and  positions  r,  the  equations 
of  motion  are  (Thee,  1986:9): 


N  _  dvi 

Y.  F  ■  =  <n  -  (  13  ) 

i  dt  , 

where  F*  is  the  force  of  atom  i  with  mass  m  and  acceleration 
dv  i  /dt .  And 


£ 


F* 


du( r  ,  ,  ) 


dr  i  j 


(  14  ) 


where  du/dr  is  the  derivative  with  respect  to  the  atom 
positions,  r,  of  the  interatomic  potential,  u.  Thus 


dv  i 
dt 


N 

£ 


j*i 


du(ru) 


dr  i  j 


(  15  ) 
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The  above  equations  were  numerically  integrated  and 
used  along  with  the  interatomic  potential  to  update  the 
atoms  positions  and  velocities.  These  positions  and 
velocities  of  the  atoms  were  updated  in  the  molecular 
dynamics  portion  of  the  code  for  a  certain  number  of  time 
steps.  These  time  steps  will  be  referred  to  throughout  the 
rest  of  this  paper  as  N3t«ps. 

The  general  format  of  the  molecular  dynamics  portion 
of  the  code  which  moves  the  atoms  Nst^ps  is  presented  below. 

R,  =  R  i  -  i  +  (  2v  ,  +  Fi  dt)(dt/2).  (16) 

The  force,  F,  is  updated  by  evaluating  the  positions 
using  an  interatomic  potential.  This  new  force,  FI,  is  used 
to  calculate  the  new  velocities.  These  new  velocites  are 
then  used  in  the  next  cycle  to  calculate  the  new  positions 
with  Eq.  (16).  The  new  velocities  are,  as  the  mass  is  set 
to  1  : 

v ,  =  v;.,  +  (F  +  FI ) ( dt/2 )  ( 17 ) 

The  Eqs .  (16)  and  (17)  were  thus  utilized  Nsteps  times 

to  move  the  atoms.  This  methodology  is  used  until  the 
atoms  have  moved  a  given  number  of  steps.  Each  step  is 
constrained  to  move  within  a  certain  period  of  time,  dt . 

Thus  the  atoms  move  N,  t«pi  with  each  step  moving  in  the 
specified  time  interval,  dt.  The  particular  outline  above 
was  used  in  the  molecular  dynamics  portion  of  the  new 
annealing  technique  to  move  the  atoms  a  prescribed  distance 
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after  a  defintive  amount  of  energy  was  applied  to  the 
system . 

Molecular  Statics. 

The  basic  purpose  of  a  molecular  statics  program  is  to 
determine  the  atomic  positions  of  a  system  such  that  the 
total  potential  energy  is  minimized.  That  is  the  minimized 
potential  energy  system  represents  stable  and  metastable 
configurations  of  the  atoms.  The  calculated  molecular 
statics  stable  and  metastable  configuration  results  are 
strictly  valid  only  at  the  temperature  of  absolute  zero 
(Sabochick,  1985:64). 

Molecular  statics  has  been  historically  divided  into 
two  general  categories.  The  first  of  these,  the  quasi- 
dynamical  technique,  is  a  method  whereby  the  Newtonian 
equations  of  motion  are  solved.  Utilizing  a  gradual  removal 
of  energy  from  the  system  and  the  continuing  solutions  of 
the  Newton  equations  over  numerous  timesteps  the  minimum 
energy  configuration  is  calculated.  The  second  type  of 
minimization  techniques  do  not  involve  directly  the 
solution  of  Newton’s  equations  of  motion.  Instead  these 
minimization  techniques  use  mathematical  methods  to  obtain 
the  minimum  energy  configurations  of  the  system.  These 
mathematical  methods  are  known  as  variational  methods 
(Sabochick,  1985:68-70). 

|  Two  common  quasidynamic  methods  are  the  frictional  and 

) 
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v.a.  methods.  The  fictional  method  allows  a  braking  force 
to  be  applied  to  the  atoms  which  is  proportional  to  the 
atoms  velocities.  The  rationale  for  this  approach  stems 
from  the  total  energy,  which  is  the  sum  of  the  potential  and 
kinetic  energies,  remaining  constant.  As  the  kinetic 
energy,  mv*/2,  increases  the  kinetic  energy  goes  to  its 
maximum  value  and  the  potential  energy  goes  to  its  minimum 
value.  Thus  the  greatest  applied  braking  force  happens  when 
the  potential  energy  is  at  a  minimum  which  slows  down  the 
system  near  the  minimum  configuration  ( Sabochick , 1 986 ) . 

The  second  method  could  be  termed  the  velocity  and 
acceleration  method  or  the  v.a  method.  This  involves 
setting  the  velocity  to  zero  whenever  the  dot  product  of  the 
velocity  and  acceleration  becomes  negative.  The  rationale 
behind  this  method  can  be  illustrated  by  considering  a  bead 
allowed  to  slide  on  a  string.  (See  Fig  5.) 


Figure  5.  Bead  on  a  Wire.  Where  point  A  is  a 
local  energy  minimum  and  point  B  is  the  global 
energy  minimum. 
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If  there  are  no  frictional  forces  the  bead  will  cycle 
along  the  above  path  indefinitely.  This  is  precisely  the 
results  of  a  molecular  dynamics  code  if  no  frictional  forces 
are  allowed.  Note  as  the  bead  passes  through  the  minimum  its 
velocity  and  acceleration  are  in  opposite  directions.  At 
this  point  setting  the  velocity  equal  to  zero  would  force 
the  bead  to  return  to  the  minimum  energy  configuration  which 
the  bead  overshot  slightly.  One  can  find  the  minimum 
configuration  by  continuing  this  process  until  for  example, 
the  maximum  force  on  any  atom  is  below  a  given  value 
(Sabochick,  1986). 

A  number  of  mathematical  methods  have  also  been  used  as 
minimization  techniques.  One  of  these,  the  method  of 
steepest  descent,  can  be  illustrated  by  using  the  bead 
example  above.  This  method  looks  at  the  position  of  the 
bead,  calculates  the  first  derivative  and  moves  the  bead  in 
the  direction  of  steepest  descent.  This  is  done  in  very 
small  incremental  steps  until  the  bead  reaches  the  minimum 
energy  configuration.  Although  accurate  this  method  is 
computationally  very  time  consuming  as  the  time  step  size 
must  be  chosen  to  be  very  small  to  prevent  the  overshooting 
of  the  minimum. 

Another  mathematical  method  is  the  Modified  Fletcher 
Powell,  MFP,  method.  This  method  is  similar  to  the  steepest 
descent  method.  The  first  derivative  is  calculated  and  the 
bead,  using  the  same  analogy  as  before,  selects  the 
direction  it  will  move  in.  However  the  bead  now  deals  with 
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a  three-dimensional  enviroment  instead  of  the  prior  two- 
dimensional  world.  Where  formerly  the  bead  saw  only  lines 
going  up  and  down  it  now  sees  hills  and  valleys.  The  MFP 
method  estimates  the  second  derivative  of  the  system.  This 
additional  information  allows  the  bead  to  make  a  two 
dimensional  selection  of  the  direction  it  should  move  in. 
This  method  is  described  in  greater  detail  by  Sabochick 
(Sabochick,  1985:72-79).  The  Modified  Fletcher  Powell 
method  was  used  as  the  statics  portion  of  the  code  to  verify 
Dagens  potentials  and  in  the  new  annealing  technique. 

Note  that  all  of  the  Molecular  Statics  techniques 
presented,  even  the  Modified  Fletcher  Powell  Method,  will 
force  the  bead  into  position  A.  Position  A  is  a  local 
energy  minimum  while  point  B  is  the  global  energy  minimum. 
The  global  energy  minimum  is  most  desirable  to  obtain  as  it 
is  the  most  stable  and  thus  the  most  likely  configuration  of 
the  system.  Finding  this  global  minimum  will  be  addressed 
in  the  annealing  sections  of  this  chapter. 

Langevin  Dynamics. 

Langevin  dynamics  was  used  by  Biswas  and  Haaman 
to  find  global  minimums  of  three-dimensional  clusters  of 
silicon  atoms.  The  utilization  of  the  Langevin  equations  is 
associated  with  the  Brownian  type  of  motion  that  the  parti¬ 
cles  experience  as  they  interact  with  each  other.  The 
equations  are  derived  from  the  assumption  of  a  very  small 
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time  variation  associated  with  a  large  fluctuating  force. 

The  motion  of  a  particle  traveling  through  the  atmosphere  or 
through  a  heated  material  has  these  characteristics  as  it 
collides  frequently  with  other  particles.  Also  at  each 
collision  a  large  force  is  transferred  to  the  particle.  As 
this  force  tends  to  be  in  the  opposite  direction  of  the 
general  forward  path  of  the  particle,  the  particle  is  slowed 
down  over  time.  This  viscous  motion  of  the  particle  with 
the  widely  varying  applied  force  is  included  in  the  Langevin 
equation  in  the  following  manner: 

dv  , 

-  =-g’  v  ,  +  F ,  ( t )  (18) 

dt 

where  g’  is  the  viscous  damping  force,  v  is  the  velocity, 

F , ( t )  is  the  time  dependent  force,  and  dv * /dt  is  the  time 
derivative  of  the  velocity  commonly  known  as  the  accelera¬ 
tion.  The  force  F,  ( t)  is  a  time  averaged  force  that 
depends  on  the  bulk  properties  of  the  system  (Nichelson, 

1983 : 267-269 ) . 

Biswas  and  Haaman  in  their  article  note  that  both  the 
equipartition  of  energy  theorem  and  the  Maxwellian  velocity 
distribution  can  be  derived  from  the  integrated  form  of  Eq . 
(16).  This  implies  that  the  particles  described  by  the 
Langevin  equations  will  achieve  an  equilibrium  temperature 
distribution.  The  BH  method  using  the  Langevin  equations  do 
not  provide  accurate  instataneous  positions  and  velocities 
of  the  particles.  Only  at  temperatures  of  absolute  zero  are 
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their  results  valid.  Also  note  that  the  Langevin  equations 
are  singly  dependent  on  temperature  in  BH ’ s  work  (BH.1986). 

Biswas  and  Haaman  modified  the  above  equations  to 
include  a  position-dependent  internal  force  h(x)  which 
yields  the  following  coupled  equations  where  the  i  indices 
have  been  dropped  for  convenience.  The  equations  are 
(Biswas  and  Haaman,  1986): 

dx 

-  =  v  (19) 

dt 

dv 

-  =  h ( x )  -  g'  v  +  F(t)  (20) 

dt 

Biswas  and  Haaman  noted  that  numerically  integrating 
the  above  differential  equations  yields  (Biswas  and 
Haaman ,  1986 ) : 

x(dt)  =  x0  +  (dt/2)  [(1  -  g’  dt)  v0  +  dt  h(y0) 

+  (dt  e)  ‘ *  Z]  (21  ) 


and 


,  v(dt)  =  v0  +  (dt/2)  [  (1  -  g’dt)  h ( y o ) 

I 

+  h(y0  +  dt  v0)  -  2g’  [  1  -  (g’  dt/2)]  v0 

-  g’  (dt  e) ‘  '»  Z  ]  +  (dt  e)  1  *  Z  (22) 

where  dt  is  the  time  step  size,  g’  is  the  viscous  damping 

f 

I 
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force,  x0  and  v0  are  respectively  the  initial  positions  and 
initial  velocities.  And  h(y0)  is  the  internal  force  depend¬ 
ing  upon  the  positions  of  the  atoms.  This  is  the  force  calcul 
ated  from  the  Stillinger  and  Weber  silicon  potential.  The  e 
is  the  product  of  2  g’kb  T,  where  kb  is  the  Boltzman 
constant  and  T  is  the  averaged  temperature  in  degrees  kelvin. 
Finally  Z  is  a  single  random  variable  derived  from  a  gausian 
with  mean  of  zero  and  standard  deviation  of  1.0  (BH.1983). 

Note  that  setting  the  viscous  damping  force,  g’,  equal 
to  zero  reduces  the  Langevin  equations  to  the  more  familiar 
form  of  molecular  dynamics  equations  presented  earlier.  The 
resetting  of  g*  was  used  to  toggle  the  molecular  dynamics 
code  from  the  Langevin  equations  to  the  earlier  set  of 
molecular  dynamics  equations. 

The  Two  Annealing  Techniques. 

Annealing  of  crystals  has  been  done  practically  in 
industry  as  a  means  of  strengthening  materials.  The 
material  is  heated  to  a  high  temperature  at  which  the 
crystalline  atoms  are  free  to  arrange  themselves  in  new 
patterns.  These  new  atomic  arrangements  tend  to  be  in  the 
patterns  prescribed  by  the  lowest  energy  minimum  the 
perfect  crystal  structure.  To  achieve  an  amorphous,  non- 
crystalline  arrangement  farthest  from  the  perfect  crystal 
structure,  the  material  is  cooled  very  rapidly  through  a 
quenching  process.  If  a  more  perfect  crystalline  structure 
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needs  to  be  achieved,  as  in  the  case  of  lenses  for 
telescopes,  another  process  is  used.  The  material  is 
softened  through  the  application  of  heat  and  then  allowed  to 
cool  in  gradual  stages.  The  gradual  nature  of  this  cooling 
process  allows  the  atoms  time  to  arrange  themselves  into 
their  preferred  states.  This  rearrangement  tends  to  form  a 
perfect  crystal  structure  (Chalmers,  1982:109-111). 

The  perfect  crystaline  atomic  arrangement  through  the 
latter  annealing  process  can  be  modeled  via  computer 
simulation.  Biswas  and  Hamaan  sugested  a  technique  for 
modeling  the  annealing  of  sililcon  utilizing  Langevin 
dynamics.  This  process  is  computational  time  intensive  as 
the  atomic  positions  must  be  calculated  200-300  times  at 
each  temperature  increment.  Sabochick  proposed  using  a 
combination  of  molecular  dynamics  and  molecular  statics  to 
model  the  annealing  process.  The  two  techniques  are 
described  below  along  with  a  number  of  computational  tools 
used  in  each  technique.  The  last  section  of  this  chapter 
will  deal  with  the  common  computational  tools  used  such  as 
the  constraining  method  and  related  parameters  used  to  hold 
the  atoms  together  (Biswas  and  Haaman ,  1986). 

Biswas  and  Haaman  Annealing  Method. 

Biswas  and  Haaman ’ s  annealing  technique  consists  of 
heating  the  crystal  to  a  melt  state  and  allowing  the  atoms 
to  reposition  themselves  at  this  plateau  temperature 
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utilizing  Langevin  molecular  dynamics. 


Then  they 


incrementally  lower  the  plateau  temperature  until  the 
system  attains  the  temperature  of  absolute  zero.  At  each 
plateau  they  allow  the  atoms  to  reposition  themselves  by 
utilizing  Langevin  dynamics.  This  gives  a  final  annealed 
configuration  that  they  believe  is  closer  to  the  global 
minimum  configuration  than  that  which  would  be  obtained 
using  molecular  statics  alone.  One  time  through  each 
temperature  plateau  calculations  will  be  refered  to  as  a 
cycle  in  this  paper. 

The  quality  of  the  annealed  configuration  depends  on 
the  number  of  incremental  temperature  plateaus  allowed.  The 
smaller  the  temperature  differential  is  between  the 
plateaus  the  more  perfect  the  crystal  configuration  formed. 
The  amount  of  time  an  atom  stays  at  the  higher  temperatures 
is  inversely  proportional  to  the  difference  in  the 
temperatures  at  each  plateau.  At  the  higher  temperatures  the 
atoms  have  enough  energy  to  avoid  being  trapped  in  higher 
energy  metastable  states.  Unfortunately,  as  the  temperature 
differential  decreases  in  size,  the  number  of  heating  and 
cooling  cycles,  Ncy,  increase.  The  differential  temperature 
or  energy  change,  Dea,  was  set  before  the  heating  and 
cooling  cycles  were  implemented.  The  relationship  was: 

Dea  =  Ea  /  Ncy  ( 23 ) 

where  Ea  is  the  plateau  temperature  distribu'  ;d  by  adding 
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Ea  amount  of  energy  to  each  atom.  At  the  end  of  each 
heating  and  cooling  cycle  the  temperature  is  decreased  per 
the  following  equation: 

Ean*w  =  Ea0id  -  Dea.  (24) 

Thus  the  greater  the  number  of  cycles  the  less  the 
temperature  differential  between  the  cycles  as  mentioned 
briefly  before.  By  increasing  the  number  of  cycles  the 
annealing  of  the  crystal  will  approach  the  global  minimum 
configuration.  Biswas  and  Haaman  in  their  article  also 
explain  that  the  random  walk  nature  of  the  Langevin 
equations  should  allow  the  atoms  to  anneal  out  of  metastable 
states.  This  results  from  the  probability  that  the  walk 
will  be  uphill  as  well  as  downhill  in  the  potential 
(Biswas  and  Hamaan ,  1986:896). 

Parameters 

BH  chose  their  time  step  size  to  be  small  such  that  g’ 

*  dt  =  0.2  as  this  value  yielded  stable  trajectories  of  the 
particles.  Their  initial  starting  temperatures  were  set  at 
0.5  eV  and  0.24  eV  etc.  in  reduced  units.  Temperatures  in 
the  range  of  approximately  1.0  -  2.0  eV  they  found  caused 
the  crystal  to  disassociate  into  individual  sets  of  small 
atom  groups  (BH.1986).  When  evaluating  the  parameters  to  be 
used  in  their  technique  with  our  code  we  found  that  Ea  of 
0.5  eV,  with  Nsteps  of  300  and  dt  of  0.07  caused  the  cluster 
to  disassociate  even  using  the  constraining  potential 
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delineated  later.  However  at  Ea  of  0.25  eV,  Nsteps  of  300, 
and  dt  of  0.1  the  cluster  of  atoms  held  together.  Thus  the 
latter  set  of  parameters  were  used  in  the  production  runs. 

Note  the  above  parameters  were  given  in  reduced  units 
where  dt  of  1  unit  is  equivalent  to  1.85  picoseconds  ( BH , 

1986:897).  Once  these  parameters  were  set  the  only  variable 
that  was  allowed  to  change  was  the  number  of  cycles  (Ncy) 
through  wh ich  the  atoms  were  allowed  to  be  heated  and 

i 

cooled.  j 

BH  Computational  Tools 

Biswas  and  Haaman  state  that  the  general  conclusions 
drawn  from  their  article  can  be  utilized  no  matter  what 
potential  is  used  to  calculate  the  interatomic  forces. 

Although  Biswas  and  Haaman  used  Stil linger  and  Weber's 
potential,  their  article  implies  that  they  modified  the 
three  body  term  by  writing  it  as  the  sum  of  a  number  of  two 
body  terms  (Biswas  and  Haaman,  1986:897).  This  modification 

•r, 

allowed  them  to  calculate  the  forces  and  energies  in  N2 
steps  instead  of  the  usual  N3  steps  required  for  the  three 
body  term.  We  used  the  three  body  Stillinger  and  Weber  (  SW  ) 
potential  with  certain  parameter  modifications  specified  in 
Chapter  Two.  The  SW  interatomic  potential  was  used  in  our 
code  for  both  techniques.  However  we  also  included  other 
modifications  to  the  code  that  reduced  the  number  of  time 
steps  required.  These  modifications  were  used  in  both  the 
new  and  BH  annealing  techniques  and  are  described  in  more 
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detail  by  Thee  (Thee,  1986:28-36). 


New  Annealing  Method. 

The  new  annealing  technique  proposed  in  the  present 
work  uses  a  combination  of  molecular  statics  and  molecular 
dynamics.  The  crystal  is  heated  so  that  the  system  is  in  a 
liquid  amorphous  state.  This  is  done  by  imparting  to  each 
atom  a  high  energy.  This  energy  impartation  to  each  of  the 
individual  atoms  will  be  referred  to,  for  convenience,  as  a 
temperature  change.  The  atoms  are  then  allowed  to  move 
according  to  the  dictates  of  the  interatomic  potential  via 
the  molecular  dynamics  code.  Then  to  bring  the  atoms  to  a 
metastable  configuration  the  statics  code  is  employed.  Note 
that  the  statics  portion  of  the  code  .also  employs  the 
interatomic  potential  to  bring  the  atoms  down  to  their 
lowest  energy  configuration.  The  temperature  is  then 
incrementally  reduced.  The  above  process  of  utilizing 
molecular  dynamics  to  move  the  atoms  and  molecular  statics 
to  reposition  the  atoms  to  their  metastable  state  is 
repeated  for  this  new  temperature.  This  cyclic  process  is 
repeated  until  the  decremented  temperature  reaches  absolute 
zero . 

The  molecular  dynamics  portion  of  the  new  method 
allows  the  atoms  the  freedom  of  motion  necessary  to  achieve 
an  annealed  crystal.  This  freedom  allows  the  atoms  to 
rearrange  themselves  in  such  a  way  as  to  avoid  being  trapped 
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in  local  minimum  configurations.  BH  essentially  used  this 
same  underlying  philosophy  in  their  technique.  However  to 
achieve  the  minimum  energy  configuration  BH  had  to  choose 
their  step  size  to  be  extremely  small  otherwise  their 
integrator  would  become  unstable  ( BH , 1 986 : 897 ) .  The  smaller 
the  step  size,  the  larger  amount  of  computer  time  required. 
As  the  statics  routine  forces  the  atoms  to  a  minimum  config¬ 
uration  in  a  very  short  period  of  time  this  new  annealing 
method  should  be  computationally  faster  than  BH’s  method. 

Parameters 

Nstepsi  Ei  and  dt  were  the  parameters  that  required 
setting  for  the  new  annealing  technique.  If  dt  is  extremely 
small  the  atoms  do  not  move  significantly  and  the  new 
technique  requires  a  large  number  of  N«iep«  to  move  the 
atoms.  Recall  that  the  MD  portion  of  the  code  moved  the 
atoms  in  such  a  way  that  when  the  statics  technique  was 
applied  a  new  minimum  configuration  was  probable  at  each  of 
the  higher  temperature  plateaus.  Thus  with  too  small  of  a 
time  interval,  dt ,  the  system  would  not  anneal  properly.  On 
the  other  hand,  if  dt  is  chosen  to  be  too  large,  the  atoms 
will  move  too  far.  This  large  movement  may  result  in 
atoms  passing  through  each  other.  As  this  process  is 
unrealistic  the  timestep  must  be  selected  such  that  the 
atoms  move  significantly  enough  to  anneal  the  crystal  while 
still  representing  the  physical  boundaries  placed  on  the 
crystal  by  classical  physics. 


34 


These  parameters,  Nst*ps, 


dt  and  E»  are  related  by  a 


simple  formula  given  below: 


dist  =  dt  N  *  t  «  p  s  E  »  1  '  2 


(25) 


where  dist  represents  the  overall  distance  the  atoms  in  the 
system  have  moved.  To  optimize  this  distance  moved  without 
the  atoms  passing  through  each  other,  the  distance  was  set 
to  1  unit  which  is  slightly  less  than  the  bondlength  of  1.12 
units.  By  maximizing  E»,  minimizing  Nsteps  and  using  the 
above  formula  to  calculate  dt  the  parameters  were  chosen. 
Thus  the  parameters  chosen  for  the  new  annealing  technique 
were  Nstepa  of  20,  Ea  of  0.5  units,  and  dt  of  0.07  units  for 
the  first  statics  runs.  Subsequent  runs  used  Nst»ps  of  20, 

Ea  of  0.25  and  dt  of  0.1  (Table  VII).  This  latter  set  of 


TABLE  VII 

Annealing  Comparison  Run  Parameters.  New 
Techl  and  Langevinl  are  respective?  runs  made  with  the 
new  technique  and  B&H ’ s  Langevin  technique  starting  at 
a  metastable  state  with  Potential  energy  of  -44.5934. 
And  similarly  the  techniques  with  2’s  represent  runs 
made  starting  at  the  metastable  state  having  a 
potential  energy  of  -43.0642. 


Parameters 

New  Techl 

Langevinl 

New  Tech2 

Langevin2  ’ 

Ea 

0.50 

0.25 

0.25 

0.25 

N  a  t  c  p  a 

20 

300 

20 

300 

Dt 

0.07 

0.10 

0. 10 

0.10 

Gamma 

0.00 

0.02 

0.00 

0.02 

Rsw 

6.21 

6.21 

7.96 

6.21 

parameters  were  used  in  the  BH  computation.  The  second  set 
of  data  was  used  in  the  new  annealing  technique  as  the 
results  for  Ea  of  0.25  was  significantly  better  than  runs 
with  Ei  of  0.5.  Note  the  above  parameters  are  given  in 
reduced  units.  Once  these  parameters  were  set  the  only 
variable  that  was  allowed  to  change  was  the  number  of 
cycles,  Ncy,  through  which  the  atoms  were  allowed  to  be 
heated  and  cooled. 

A  Procedure  Used  in  Both  Annealing  Methods. 

Constraining  Potential. 

We  found  that  the  addition  of  a  constraining  potential 
subroutine  was  necessary  in  the  molecular  dynamics  portion 
of  both  techniques.  Without  the  addition  of  this  restraint 
the  results  of  the  BH  Langevin  and  new  annealing  method’s 
atomic  configurations  disasociated .  That  is  the  resulting 
system  had  fragmented  into  smaller  clusters  of  atoms  when  an 
Ei  of  0.5  or  0.25  eV  was  applied.  The  constraining 
potential  was  derived  such  that  within  a  certain  fixed 
radius  the  constraining  force,  FCoa,  and  constraining 
energy,  Ecoa,  terms  were  negligible.  Outside  the  chosen 
radius  Ecoa  and  FCoa  were  large  enough  to  effectively  return 
any  straying  atoms  back  within  the  circular  boundary  defined 
by  the  radius.  This  radius  was  designated  as  the  switch 
radius,  Rsw,  and  was  used  as  a  toggle  parameter  distinguish- 
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ing  between  the  negligible  and  large  constraining  potentials, 
The  energy  and  force  equations  are: 


=  Li  r 2 


r  <  rsw  ( 26  ) 


Econ  =  ao  +  a i  +  Li  r! 


r  >  rsw  ( 27  ) 


2L  i  ( X  i  -  x  a 


r  <  rsw 


Fcob  =  -  ( a i /r  +  2L 2 ) (x,  -  xiv4) 


+  (l/N)^(ai/r)  (Xi  -  x  a  v  e  )  r  >  rsw  (29) 


where  Ec<i  and  FCon  are  respectively  the  constraint  energy 
and  force.  Li  and  L2  are  respectively  the  inner  and  outer 
radius  constraint  parameters.  The  assigned  values  of  these 
two  parameters  determine  the  strength  of  the  force  field 
applied  to  the  atoms  under  their  jurisdiction.  As  we 
require  a  negligible  force  within  Rsw,  Li  is  equal  to 
0.0001.  And  L2  is  equal  to  0.1  as  a  fairly  large  force  is 
required  outside  of  the  radius  Rsw.  N  in  the  above  equation 
is  the  number  of  particles,  a0  and  ai  are  constants  selected 
that  the  toggle  radius  values  are  continuous  and  smooth.  The 
constraints  are: 


a  i  =  2 ( L  i  -  L  2 )  Rsw 


a  2  =  -  (Li  -  L  2 )  Rsw  * . 
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And  ri  is  the  radial  distance  given  by  the  following 
equation 


r,  =  [(Xi  -  x,„)!  +  (yi  -  ya»«)2]l/2  (32) 

where  x,  and  y,  are  the  two-dimensional  position  of  atom  i 
in  rectangular  coordinates.  The  center  of  the  circular 
system  is  represented  by  xav*  and  ya«{.  The  equations  are: 


Xa  v  ( 


(  33) 


Note  that  Eqs .  (28),  (29)  and  (33)  only  present  the  x 

component  of  the  force.  These  equations  would  thus  have  to 
be  repeated  for  the  y  components. 

Rsw  was  calculated  for  all  but  the  second  BH  Langevin 
runs  by  calculating  the  radius  of  the  outer  most  atom  from 
the  center  of  the  initial  configuration  used.  Then  a  buffer 
zone  of  1.0  was  added  to  this  calculated  value.  For  the 
second  Langevin  run  the  value  of  6.21  was  hand  set  to 
provide  a  tighter  boundary. 

The  above  equations  were  used  every  time  the  forces  and 
energies  were  calculated  via  the  interatomic  potentials.  If 
the  atom  was  within  the  prescribed  radius  Rsw  the  change  in 
the  total  energy  and  force  of  the  atom  cluster  was 
negligible.  If  on  the  other  hand,  the  atom  was  beyond  the 
prescribed  radius  the  interatomic  energy,  E,  and  the 
interatomic  force,  F ,  ,  terms  were  updated  with  a  non- 
negligible  additive  term. 
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The  equations  used  to  update  the  potential  energy  and 
force  are: 

E  =  E  + ^ E  c  o  o  (31) 

and 

F,  =  F,  +  Fcon.  (35) 

Note  once  again  only  the  x  component  of  the  force  is 
given.  A  y  component  equation  is  also  neccessary  to  model 
the  two-dimensional  36  atom  cluster  that  was  chosen  as  our 
test  case. 


IV.  SILVER  CALCULATION 


Dagens'  Potential  Calculations. 

The  purpose  of  the  first  part  of  this  thesis  is  to 
verify  whether  two  new  potentials  by  Dagens  are  more  accurate 
than  an  older  potential.  The  basis  for  this  comparison  is 
the  experimental  data  value  of  1.11  eV  for  the  vacancy 
energy  of  formation  for  silver  (Lam,  1983:2514).  The  value 
for  the  vacancy  energy  of  formation  was  computationally 
calculated  for  the  three  potentials  and  compared  to  the 
above  experimental  data  value. 

To  validate  Dagens’  potentials  the  Modified  Fletcher 
Powell,  MFP,  method  of  molecular  statics,  MS,  was  employed. 
This  method  was  used  as  a  working  molecular  statics  code  was 
available.  In  particular,  molecular  statics  was  used  in 
conjunction  with  Dagens’  potentials  to  calculate  the  energy 
of  formation  of  a  single  vacancy. 

The  energy  of  formation  of  a  vacancy  was  calculated  by 
first  utilizing  a  perfect  crystal  with  no  defects.  The 
minimum  energy  of  the  perfect  crystal  was  then  determined  by 
using  the  MFP  MS  code.  Next  one  or  more  of  the  inner  atoms 
were  removed,  forming  a  single  or  multiple  vacancy.  The 
minimum  energy  of  this  modified  crystal  was  then  calculated. 
These  minimum  energies  correspond  to  the  stable  or 
metastable  states  of  the  system.  To  calculate  the  energy  of 
formation  of  a  vacancy,  the  following  mathematical  equation 
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was  employed: 
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where  the  first  term  is  the  minimum  energy  of  the  perfect 
crystal  found  by  using  the  molecular  statics  technique 
previously  presented.  The  second  term  summation  is  the 
minimum  potential  energy  of  a  crystal  with  a  vacancy,  an 
inner  atom  removed.  This  summation  term  is  multiplied  by  a 
modifying  term.  In  this  term,  n  is  the  number  of  atoms  in 
the  perfect  crystal  and  v  is  the  number  of  vacancies  in  the 
modified  crystal.  The  last  term  is  the  virial  term  which 
compensates  for  the  changed  volume  of  the  system.  The 
virial  term  is  calculated  during  the  perfect  crystal  energy 
minimization  process  (Lara,  1983). 


Code  Verification. 

To  verify  that  the  code  was  working  properly  several 
calculated  values  were  compared  with  Lam’s  values  for 
various  radius  dimensions.  The  operational  computer  code 
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was  verified  by  calculating  the  vacancy  energy  of  formation 
for  the  old  potential.  This  value  was  compared  to  Lam’s 
published  value  (Lam, 1983;  Lam, 1986).  This  process  resulted 
in  close  correlation  with  the  values  Lam  obtained.  The 
results  of  this  check  is  presented  in  Table  VIII.  Lam’s 
data  values  were  based  on  a  864  atom  ensemble  with  a  cut  off 
radius  of  11.632  A0  (Lam, 1983).  The  calculated  values  were 
determined  using  a  256  atom  ensemble  with  a  cut  off  radius 
of  6.807  A0 . 


TABLE  VIII 

Multiple  Vacancy  Formation  Energies.  The  table 
compares  Lam's  Values  for  the  old  Dagen  potential 
values  with  this  projects  calculated  values. 


VACANCIES 

LAM’S  VALUES 

CALCULATED 

VALUES 

SINGLE-VACANCIES 

1  .37 

1  .  34 

DI-VACANCIES 

2 . 70 

2 . 08 

TRI-VACANCIES 

3  .65 

3  .  77 

TETRA- VACANCIES 

4.62 

4.68 

The  cut-off  radius  is  another  parameter  that  can  be 
used  as  a  rough  comparison  to  gauge  the  accuracy  of  the 
code.  The  cut-off  radius  is  the  minimum  distance  the  atoms 
can  be  from  one  another  while  still  experiencing  a 
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negligible  force  from  the  interation  (Fig  6).  Atoms  outside 
the  cut-off  radius  of  a  particular  atom  do  not  affect  that 
particular  atom.  Also  atoms  that  are  within  the  cut-off 
radius  of  a  particular  atom  are  influenced  by  that  atom  and 
themselves  affect  the  atom.  The  cut-off  radius  for  this 
paper  was  taken  to  be  a  few  lattice  constants  away  from  any 
particular  atom.  Lam  also  chose  the  cut-off  radius  for 
Dagens  silver  potentials  to  be  at  least  twice  the  silver 
lattice  constant  of  4.087  A0.  In  addressing  the 
verification  of  the  old  Lam  potential  the  cut-off  radius  was 
set  and  then  refined  during  the  calculation  of  the  vacancy 
formation  energies.  Lam  calculated  the  value  of  the  cut-off 
radius  for  the  864  atom  system  to  be  11.737  A0  (Lam, 1983). 

My  calculated  values  were  6.807  A0  for  the  256  atom  system 
and  12.248  A0  for  the  864  atom  system. 


Figure  6.  Cut  off  Radius.  Illustrated  with  a 
pair  potential  for  two  isolated  atoms. 


43 


Results  and  Discussion 

The  computational  tools  of  molecular  statics  were  used 
to  verify  Dagens ’  two  new  potentials.  The  results  obtained 
are  presented  Tables  IX  and  X. 


TABLE  IX 

Vacancy  Formation  Energy  of  the  Silver 
Potentials.  (256  atom  ensemble) 


cut  off 

Virial 

Formation 

Ref  erence 

radius  (A) 

(eV) 

Energy  ( eV ) 

Actual  Experi- 

mental  Results 

N/A 

N/A 

1.11 

Lam’s  Results 

(old  potential) 

11.737 

-2.096 

1  .  37 

Calculated  Value 

( old  potential ) 

6 . 807 

-1.991 

1.34 

Calculated  Value 

( new  potential  1 ) 

6 .923 

-0.716 

0 . 58 

Calculated  Value 

(new  potential  2) 

7.014 

-0.075 

0 . 37 

The  experimental  value  is  the  main  comparison  value. 
This  value,  1.11  eV  represents  the  energy  of  formation  of  a 
single  vacancy  in  a  perfect  crystal.  The  removed  atom  is 
assumed  to  have  moved  from  the  interior  of  the  system  to 
infinity  thus  the  256  and  864  atom  ensembles  should  yield 
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TABLE  X 


Vacancy  Formation  Energy  of  the  Silver 
Potentials.  (864  atom  ensemble; 


cut  off 

Virial 

Formation 

Reference 

radius  (A) 

( eV ) 

Energy  (eV) 

Actual  Experi- 

mental  Results 

M/A 

M/A 

1.11 

Lam’s  Results 

( old  potential ) 

1 1  .  737 

-2 . 096 

1  .37 

Calculated  Value 

( old  potential ) 

12.248 

-2. 109 

1.41 

Calculated  Value 

( new  potential  1 ) 

6.923 

-0.643 

0.51 

Calculated  Value 

(new  potential  2) 

9 . 869 

-0.039 

0.31 

approximately  the  same  results.  The  energy  of  formation 
is  one  of  the  first  parameters  commonly  checked  when 
validating  a  potential  and  thus  was  selected  as  the  first 
verification  point  for  the  Dagens 1  potentials  under 
consideration. 

The  calculated  values  in  tabular  form  on  the  preceding 
pages  are  the  values  obtained  using  the  molecular  statics 
technique  for  respectively  a  256  atom  three-dimensional 
ensemble  and  a  864  atom  three-dimensional  ensemble.  Lam’s 
values  for  his  old  potential  were  calculated  utilizing  a  864 
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atom  three-dimensional  ensemble  and  are  as  noted  above  also 


approximately  valid  for  the  256  atom  three-dimensional 
ensemble . 

The  last  column  of  data  in  each  table  is  the  comparison 
values .  Note  that  the  two  new  potential  values  for  the  256 
atom  ensembles  vary  respectively  by  48%  and  67%  from  the 
experimental  data  value.  While  my  calculated  value  of  the 
old  potential  varies  by  only  21%.  The  864  atom  ensemble 
yields  similiar  results.  Thus  the  old  Dagens  potential 
portrays  more  accurately  the  vacancy  formation  energy  than 
the  two  new  potentials  do.  As  the  vacancy  formation  energy 
for  the  new  potentials  is  less  accurate  than  the  previous 
potentials  value,  the  conclusion  was  reached  that  the 
previous  Dagens  potential  should  not  be  replac-ed  by  either  of 
the  two  new  Dagens  potentials. 
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IV.  ANNEALING  CALCULATIONS 


The  purpose  of  the  second  part  of  this  thesis  was  to 
investigate  a  computationally  faster  method  of  simulating 
the  annealing  of  atomic  clusters.  This  second  task  required 
the  comparison  of  a  technique  presented  in  an  article  by 
Biswas  and  Haaman  and  a  new  anealing  technique.  The  Biswas 
and  Haaman,  BH ,  method  utilized  Langevin  dynamics  to  anneal 
the  material  (Biswas  and  Haaman,  1986)  .  The  new  method  used 
a  combination  of  molecular  dynamics  and  molecular  statics  to 
simulate  the  annealing  of  the  material  in  question.  These 
two  techniques  were  compared  on  the  basis  of  which  method 
achieved  a  global  minimum  or  near  metastable  state 
computationally  faster  than  the  other  procedure.  Both 
methods  were  evaluated  using  a  modified  form  of  the 
Stillinger  and  Weber,  SW ,  silicon  potential.  In  general  the 
atoms  were  thermally  excited  to  the  melt  point  and  then 
allowed  to  cool  through  a  series  of  heating/  cooling  stages 
to  a  global  minimum  configuration  (Biswas  and  Haaman,  1986), 
Our  goal  in  comparing  the  two  techniques  was  to  find 
the  method  which  reached  a  particular  minimum  in  the  fastest 
time.  Also  required  was  that  the  atoms  anneal  to  their 
global  minimum  configuration  or  a  metastable  configuration 
near  the  global  minimum.  The  major  portion  of  the  code  that 
both  the  statics  and  two  dynamics  routines  frequently  use  is 
the  calculation  of  the  interatomic  forces  and  energies  as  a 
function  of  the  atoms’  positions.  This  subroutine,  Eval, 
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required  a  large  portion  of  computer  time  as  it  calculated 
the  energies  and  forces  using  the  modified  SW  potential. 
Keeping  track  of  the  number  of  times  this  subroutine  is 
called  will  give  a  fair  estimate  of  the  relative  speeds  of 
the  two  methods.  Also  as  the  two  methods  were  tested  on  the 
ASD  cray  computer,  actual  CPU  times  are  available  and 
provide  another  data  point  in  estimating  the  computational 
intensity  of  the  two  methods. 

The  simulation  crystal  was  chosen  to  be  a  36  atom 
cluster.  It  was  determined  to  simulate  this  crystal  in  two 
dimensions  as  a  two  dimensional  arrangement  is  easier  to 
plot  and  visualize  than  a  three-dimensional  arrangement. 

Also  a  two  dimensional  arrangement  allowed  us  to  hand 
calculate  the  minimum  configuration.  This  was  accomplished 
by  hand  drawing  the  'minimum  energy’  configuration  of  the  36 
atom  cluster.  Then  counting  the  bondlengths  and  assigning  a 
value  of  1.0  unit  to  each  of  them.  Through  this  process  the 
minimum  energy  configuration  value  of  -46.0  was  obtained 
( See  Fig  7  )  . 

Listed  below  are  the  differences  and  similarities  in 
the  computational  implemetat ion  of  both  the  BH  and  the  new 
technique.  First  the  parameters  for  each  technique  varied 
(Table  VII).  Next  for  the  Biswas  and  Haan.an  technique, 
Langevan  dynamics  were  used  in  place  of  the  usual  molecular 
dynamics  to  anneal  the  material.  Molecular  statics  was  not 
used  in  the  Biswas  and  Haaraan  technique.  Finally,  the  new 
method  employed  the  use  of  second  order  molecular  dynamics 
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and  modified  Fletcher  Powell  molecular  statics  to  anneal  the 
material.  The  code  modification  that  was  used  commonly  for 
both  annealing  techniques  was  the  constraining  system  bound¬ 
aries  designated  by  Rsw  that  contained  the  system  within  a 
circular  area.  Both  methods  were  run  using  a  two-dimensional 
crystal  of  36  atoms.  In  both  procedures  once  the  parameters 
were  set,  the  variable  used  to  achieve  a  better  anneal  was 
the  increasing  of  the  number  of  heating  and  cooling  cycles 
(Ncy).  Other  than  the  above,  both  techniques  utilized  the 
same  code,  therefore  a  direct  comparison  of  the  results  from 
the  two  procedures  can  be  made. 

Two  sets  of  runs  were  made  to  test  the  two-dimensional 
thirty-six  atom  clusters.  Each  set  of  runs  was  accomplished 
using  both  techniques  with  the  starting  positions  differing. 
The  new  annealing  runs  were  conducted  in  such  a  way  that 
after  each  temperature  decrement  a  metastable  state  of  the 
atom  configuration  was  found  through  the  MS  portion  of  the 
code.  Thus  the  minimum  energy  states  tabulated  are  the 
lowest  energy  minimum  found  by  the  MS  routine.  This  minimum 
energy  configuration  is  not  necessarily  the  final  state  of 
the  atoms  at  the  end  of  the  run.  In  the  BH  Langevin  method, 
however,  BH  specify  that  the  final  results  at  the  end  of  the 
BH  run  is  the  minimum  energy  results  obtained.  These  are 
the  BH  results  tabulated. 

The  first  set  of  runs  started  from  the  minimum  energy 
configuration  of  -44.59  eV.  The  results  of  these  runs  is 
shown  in  Table  XI  and  in  graphical  form  in  Figures  8  and  9. 
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TABLE  XI 


Set  1  Annealing  Runs  Results.  New  annealing 
Technique  and  BH  Langevin  Technique  results  starting 
from  the  minimum  energy  configuration  of  -44.593  eV . 


Ncy 

New  annealing 
Technique 

BH  Langevin 
Technique 

10 

Emin  -  -43.762 

CPU  =  110  sec 

Eval  =  6000 

20 

1 

Emin  =  -43.064 

CPU  =  193  sec 

Eval  =  12000 

40 

Emin  =  -45.3141 

CPU  =  36  sec 

Eval  =  2600 

Emin  =  -43.871 

CPU  =  359  sec 

Eval  =  24000 

80 

1 

1 

160 

1 

The  second  set  of  runs  started  from  the  minimum  energy 
configuration  of  -43.064.  The  configuration  used  was  the 
results  of  the  20  cycle  BH  Langevin  technique  in  the  first 
set  of  runs.  The  purpose  of  using  the  prior  BH  results  was 
to  answer  the  question,  'Can  the  two  methods  achieve  the 
-45.3  minimum  state  from  the  -43.064  state?’  The  results 
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Figure  8.  Set  1  Annealing  Runs  Results.  Results  starting 
from  the  minimum  energy  configuration  of  -44.593  eV. 
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of  this  new  set  of  runs  is  portrayed  in  Table  XII  and  in 
Figures  10  and  11  on  the  following  pages. 


TABLE  XII 


Set  2  Annealing  Runs  Results.  New  annealing 
Technique  and  BH  Langevin  Technique  results 
starting  from  the  minimum  energy  configuration 

of  -43.064  eV. 


Figure  10.  Set  2  Annealing  Runs  Results.  Results  starting 
from  the  minimum  energy  configuration  of  -43.064  eV . 
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STORTING  POSITION 


Figure  11.  Set  2  Annealing  Runs  Results.  Results  starting 
from  the  minimum  energy  configuration  of  -43.064  eV. 


-45.60 


Note  that  the  important  comparison  parameters  in  Tables 
XI  and  XII  are  the  number  of  evaluations  (Eval)  and  the  CRAY 
CPU  time.  The  hand  calculated  minimum  energy  configuration 
gave  a  value  of  -46.0  Both  the  new  annealing  technique  and 
the  BH  Langevin  method  achieved  -45.3  as  their  calculated 
energy  minimums  in  the  first  set  of  runs.  The  important 
point  to  note  however  is  the  amount  of  computer  time  it  took 
each  of  the  two  methods  to  achieve  this  minimum  energy 
configuration.  Where  it  took  the  new  annealing  technique  36 
seconds  of  CPU  time  and  2600  calls  to  Eval;  the  BH  Langevin 
technique  reached  this  minimum  energy  configuration 
somewhere  between  359-715  seconds  and  24000-48000  calls  to 
Eval.  Thus  the  amount  of  computer  time  required  to  reach 
the  same  local  minimum  was  approximately  a  factor  of  ten 
difference  between  the  two  techniques. 

The  next  set  of  runs  were  used  as  another  comparison 
point.  Since  the  starting  configuration  was  at  a  higher 
energy  minimum  than  before  it  was  expected  that  it  would 
take  longer  to  obtain  the  local  minimum.  From  Fig.  9  and 
the  data  in  Table  XII  one  can  see  that  the  new  technique  once 
again  achieved  the  local  minimum  faster  than  the  BH  Langevin 
method.  Due  to  a  lack  of  time  the  BH  Langevin  runs  in  set  2 
were  not  extended  to  the  minimum  energy  configuration  value. 

A  number  of  two-dimensional  plots  of  the  atom 
configurations  accompanying  the  tabulated  energy  minimum 
values,  Emin,  are  shown  on  the  following  pages.  Note  that  a 
slight  change  in  the  minimum  energy  obtained  changes  the 
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Figure  12.  Set  1,  New  Annealing  Technique,  cycle  40. 

New  annealing  Technique  results  starting  from  the 
minimum  energy  configuration  of  -44.593  eV .  The  above 
plot  corresponds  to  a  minimum  energy  configuration  of 
-45.314  and  took  36  seconds  of  CRAY  CPU  time. 


) 
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Figure  13.  Set  1,  New  Annealing  Technique,  cycle  on. 

New  annealing  Technique  results  starting  from  the 
minimum  energy  configuration  of  -44.593  oV .  The  above 
plot  corresponds  to  a  minimum  energy  configuration  of 
-45.248  and  took  69  sec  of  CRAY  CPU  time. 


Figure  14.  Set  1,  BH  Langevin  Technique,  cycle  20. 

BH  Langevin  Technique  results  starting  from  the 
minimum  energy  configuration  of  -44.593  eV .  The  above 
plot  corresponds  to  a  minimum  energy  configuration  of 
-43.064  and  took  193  sec  of  CRAY  CPU  time. 
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Figure  15.  Set  1,  BH  Langevin  technique,  cycle  80. 

BH  Langevin  Technique  results  starting  from  the 
minimum  energy  configuration  of  -44.593  eV.  The  abo 
plot  corresponds  to  a  minimum  energy  configuration 
-45.287  and  took  715  sec  of  CRAY  CPU  time. 


^20 


Figure  17.  Set  2,  BH  Langevin  Technique,  cycle  20. 

BH  Langevin  Technique  results  starting  from  the 
minimum  energy  configuration  of  -43.064  eV.  The  above 
plot  corresponds  to  a  minimum  energy  configuration  of 
-41.721  and  took  188  sec  of  CRAY  CPU  time. 
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atomic  arrangement  of  the  atoms  somewhat.  N'ote  also  from 
Fig  16  that  the  new  technique  was  able  to  take  the  higher 
energy  state  back  down  to  -45.3  in  287  seconds  of  CPU  time. 
But,  from  Fig  17.,  when  we  ran  the  BH  Langevin  technique  the 
atoms  separated  into  two  clusters.  By  lowering  the  initial 
energy  and  reaccomplishing  the  run  this  disassociation  may 
be  eliminated.  Due  to  a  lack  of  time  we  were  unable  to 
continue  these  second  set  of  runs,  but  the  above  assumption 
seems  to  be  valid  as  the  -43.064  eV  state  was  probably 
passed  through  on  the  80  cycle  BH  Langevin  run  in  Set  1 . 

From  the  small  number  of  data  points  available  the  new 
annealing  technique  seems  to  be  ten  times  faster  than  the  BH 


Langevin  method. 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  purpose  of  this  study  was  twofold:  to  determine 
the  validity  of  two  new  silver  potentials  derived  by  Dagen 
and  to  investigate  a  new  technique  for  finding  global 
minimum  configurations.  Both  of  these  projects  required  the 
use  of  the  atomistic  simulation  methods  of  molecular  statics 
and  molecular  dynamics. 

First,  molecular  statics,  MS,  was  used  to  test  two  new 
Dagens ’  silver  potentials.  In  particular  the  formation 
energy  of  a  single  vacancy  was  calculated  using  Modified 
Fletcher  Powell,  MFP,  MS  in  conjunction  with  the  new  Dagens 
potentials.  The  two  new  potential  calculated  values 
obtained  for  the  256  three-dimensional  atom  ensemble 
differed  -by  48%  and  67%  from  the  experimental  vacancy 
formation  energy  value  of  1.11  eV .  Similarly  the 
calculated  new  potential  values  obtained  for  the  864  three- 
dimensional  atom  ensemble  differed  by  54%  and  72%  from  the 
experimental  value.  The  calculated  vacancy  energy  of 
formation  value  for  the  previous  Dagens  potential  was  found 
to  vary  by  21%  for  the  256  atom  ensemble  and  27%  for  the  864 
atom  ensemble.  These  results  indicate  that  the  old 
potential  correlates  more  precisely  with  experimental  data 
than  the  two  new  potentials.  Thus  it  is  recommended  that 
the  older  version  of  the  Dagens  silver  potential  continue  to 
be  used  until  a  more  accurate  version  is  derived. 

Second  a  new  technique  for  simulating  the  annealing  of 


65 


a  material  was  investigated,  using  a  combination  of 
molecular  statics,  MS,  and  molecular  dynamics,  MD .  The 
results  of  this  technique  was  compared  with  a  prior  method 
proposed  by  Biswas  and  Haaman  utilizing  Langevin  molecular 
dynamics  techniques  (Biswas  and  Haaman,  1986).  In 
particular,  the  global  minimum  configuration  of  a  two- 
dimensional  36  atom  cluster  using  a  modified  version  of 
Stillinger  and  Weber’s  silicon  potential  was  determined. 

The  resulting  computational  times  of  the  BH  Langevin  and  the 
new  method  were  then  compared.  The  resulting  data  showed 
that  while  it  took  36  seconds  of  CRAY  CPU  time  for  the  new 
annealing  technique  to  reach  a  local  minimum  of  -45.3  eV , 
the  BH  Langevin  method  required  between  359  to  715  seconds 
of  CRAY  CPU  time  to  reach  the  same  local  minimum.  From 
these  initial  data  runs  it  appears  that  the  new  technique  is 
ten  times  faster  than  the  BH  Langevin  technique.  Thus  it  is 
recommended  that  further  investigation  be  made  into  this  new 
method  of  simulating  the  annealing  of  materials. 

The  following  itemizes  a  number  of  recommendations  for 
further  investigation  of  the  new  annealing  technique. 
Parameters  of  both  techniques  need  to  be  optimized  more 
fully.  The  BH  technique  was  implemented  in  the  last  few 
weeks  of  the  term  and  parameters  were  chosen  mainly  from 
Biswas  and  Haamans  article.  A  larger  two-dimensional  cluster 
should  be  annealed  as  this  would  indicate  what  happens  at 
grain  boundaries.  Finally  the  code  should  be  extended  to 
three  dimensions  and  the  two  techniques  compared. 
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In  conclusion,  two  new  silver  potentials  derived  by 
Dagens  were  tested  against  experimental  data  and  found  to  be 
less  precise  than  the  older  version  of  the  Dagen  silver 
potential.  Also  a  new  annealing  technique  was  compared 
with,  and  found  to  be  ten  times  faster  than  the  BH  Langevin 
annealing  method.  Investigation  of  this  new  annealing 
method  should  be  persued. 
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